
namespace g2 {

	namespace {

		// Mysterious threshold on sine of the angle between lines/planes below which they are
		// considered colinear/coplanar
		const double ALPHA = 1e-3;

	} // namespace

/*
	This code is ported from Align's Boolean2 library; 
	May contain errors, because not all functions are used and tested so far
*/

template <typename T>
bool project_point_on_segment(const basic_vector<T> & p, 
	const basic_vector<T> & p0, const basic_vector<T> & p1, 
	bool & degenerate, T & h, T & t)
{
	basic_vector<T> a = p - p0, b = p1 - p0;
	// Lengths of the segment
	T b2 = square(b);
	T eps = std::numeric_limits<T>::epsilon();
	if (degenerate = (b2 < eps * eps))
		return false;
	// Parameter of the projection
	t = (a * b) / b2;
	// Distance to the projection point
	h = (a ^ b) / sqrt(b2);
	return t >= 0 && t <= 1;
}


template <typename T>
bool intersect_segments(const basic_vector<T> & p0, const basic_vector<T> & p1,
	const basic_vector<T> & q0, const basic_vector<T> & q1,
	bool & collinear, T & s, T & t)
{
	/*
		Intersection point p must satisfy the following equations:
			p = p0 + s (p1 - p0),
			p = q0 + t (q1 - q0).
		We rewrite this as follows:
			pp.x * s - qq.x * t = - pq.x,
			pp.y * s - qq.y * t = - pq.y,
		where pp = (p1 - p0), qq = (q1 - q0), pq = (p0 - q0).
	*/
	basic_vector<T> pp = p1 - p0, qq = q1 - q0, pq = p0 - q0;

	T det = - pp.x * qq.y + pp.y * qq.x;
	
	// Changed from 16 * EPSILON to ALPHA; guessing again!
	collinear = (fabs(det) < ALPHA * abs(pp) * abs(qq)); 

	if (collinear)
	{
		basic_vector<T> n;

		n = basic_vector<T>(- pp.y, pp.x);
		n.makeUnit();

		double dq0 = (q0 - p0) * n, dq1 = (q1 - p0) * n;
		if (dq0 < -EPSILON && dq1 < -EPSILON 
			|| dq0 > EPSILON && dq1 > EPSILON)
			return false;

		n = basic_vector<T>(- qq.y, qq.x);
		n.makeUnit();

		double dp0 = (p0 - q0) * n, dp1 = (p1 - q0) * n;
		if (dp0 < -EPSILON && dp1 < -EPSILON 
			|| dp0 > EPSILON && dp1 > EPSILON)
			return false;

		Vector2D L = pp;
		L.makeUnit();
		dp0 = L * p0;
		dp1 = L * p1;
		dq0 = L * q0;
		dq1 = L * q1;
		if (dp0 > dp1)
			std::swap(dp0, dp1);
		if (dq0 > dq1)
			std::swap(dq0, dq1);
		return (dp0 < dq1 && dq0 < dp1);
	}

	s = (pq.x * qq.y - pq.y * qq.x) / det;
	t = (- pp.x * pq.y + pp.y * pq.x) / det;

	return (s >= 0 && s <= 1) && (t >= 0 && t <= 1);
}

} // g2
